////////////////////////////////////////////////////////////////////////////////////////
//
//  Copyright 2019 Alexander Stukowski
//
//  This file is part of OVITO (Open Visualization Tool).
//
//  OVITO is free software; you can redistribute it and/or modify it either under the
//  terms of the GNU General Public License version 3 as published by the Free Software
//  Foundation (the "GPL") or, at your option, under the terms of the MIT License.
//  If you do not alter this notice, a recipient may use your version of this
//  file under either the GPL or the MIT License.
//
//  You should have received a copy of the GPL along with this program in a
//  file LICENSE.GPL.txt.  You should have received a copy of the MIT License along
//  with this program in a file LICENSE.MIT.txt
//
//  This software is distributed on an "AS IS" basis, WITHOUT WARRANTY OF ANY KIND,
//  either express or implied. See the GPL or the MIT License for the specific language
//  governing rights and limitations.
//
////////////////////////////////////////////////////////////////////////////////////////

#include <ovito/grid/Grid.h>
#include <ovito/grid/objects/VoxelGrid.h>
#include <ovito/mesh/surface/SurfaceMesh.h>
#include <ovito/stdobj/simcell/SimulationCellObject.h>
#include <ovito/stdobj/properties/PropertyAccess.h>
#include <ovito/core/dataset/DataSet.h>
#include <ovito/core/app/Application.h>
#include <ovito/core/dataset/pipeline/ModifierApplication.h>
#include "CreateIsosurfaceModifier.h"
#include "MarchingCubes.h"

namespace Ovito { namespace Grid {

IMPLEMENT_OVITO_CLASS(CreateIsosurfaceModifier);
DEFINE_PROPERTY_FIELD(CreateIsosurfaceModifier, subject);
DEFINE_PROPERTY_FIELD(CreateIsosurfaceModifier, sourceProperty);
DEFINE_REFERENCE_FIELD(CreateIsosurfaceModifier, isolevelController);
DEFINE_REFERENCE_FIELD(CreateIsosurfaceModifier, surfaceMeshVis);
SET_PROPERTY_FIELD_LABEL(CreateIsosurfaceModifier, sourceProperty, "Source property");
SET_PROPERTY_FIELD_LABEL(CreateIsosurfaceModifier, isolevelController, "Isolevel");

/******************************************************************************
* Constructs the modifier object.
******************************************************************************/
CreateIsosurfaceModifier::CreateIsosurfaceModifier(DataSet* dataset) : AsynchronousModifier(dataset)
{
	setIsolevelController(ControllerManager::createFloatController(dataset));

	// Create the vis element for rendering the surface generated by the modifier.
	setSurfaceMeshVis(new SurfaceMeshVis(dataset));
	surfaceMeshVis()->setShowCap(false);
	surfaceMeshVis()->setSmoothShading(true);
	surfaceMeshVis()->setObjectTitle(tr("Isosurface"));
}

/******************************************************************************
* Determines the time interval over which a computed pipeline state will remain valid.
******************************************************************************/
TimeInterval CreateIsosurfaceModifier::validityInterval(const PipelineEvaluationRequest& request, const ModifierApplication* modApp) const
{
	TimeInterval iv = AsynchronousModifier::validityInterval(request, modApp);
	if(isolevelController()) 
		iv.intersect(isolevelController()->validityInterval(request.time()));
	return iv;
}

/******************************************************************************
* Asks the modifier whether it can be applied to the given input data.
******************************************************************************/
bool CreateIsosurfaceModifier::OOMetaClass::isApplicableTo(const DataCollection& input) const
{
	return input.containsObject<VoxelGrid>();
}

/******************************************************************************
* This method is called by the system when the modifier has been inserted
* into a pipeline.
******************************************************************************/
void CreateIsosurfaceModifier::initializeModifier(ModifierApplication* modApp)
{
	AsynchronousModifier::initializeModifier(modApp);

	// Use the first available voxel grid from the input state as data source when the modifier is newly created.
	if(sourceProperty().isNull() && subject().dataPath().isEmpty() && Application::instance()->executionContext() == Application::ExecutionContext::Interactive) {
		const PipelineFlowState& input = modApp->evaluateInputSynchronous(dataset()->animationSettings()->time());
		if(const VoxelGrid* grid = input.getObject<VoxelGrid>()) {
			setSubject(PropertyContainerReference(&grid->getOOMetaClass(), grid->identifier()));
		}
	}

	// Use the first available property from the input grid as data source when the modifier is newly created.
	if(sourceProperty().isNull() && subject() && Application::instance()->executionContext() == Application::ExecutionContext::Interactive) {
		const PipelineFlowState& input = modApp->evaluateInputSynchronous(dataset()->animationSettings()->time());
		if(const VoxelGrid* grid = dynamic_object_cast<VoxelGrid>(input.getLeafObject(subject()))) {
			for(const PropertyObject* property : grid->properties()) {
				setSourceProperty(VoxelPropertyReference(property, (property->componentCount() > 1) ? 0 : -1));
				break;
			}
		}
	}
}

/******************************************************************************
* Creates and initializes a computation engine that will compute the
* modifier's results.
******************************************************************************/
Future<AsynchronousModifier::EnginePtr> CreateIsosurfaceModifier::createEngine(const PipelineEvaluationRequest& request, ModifierApplication* modApp, const PipelineFlowState& input)
{
	if(!subject())
		throwException(tr("No input voxel grid set."));
	if(subject().dataClass() != &VoxelGrid::OOClass())
		throwException(tr("Selected modifier input is not a voxel data grid."));
	if(sourceProperty().isNull())
		throwException(tr("Please select an input field quantity for the isosurface calculation."));

	// Check if the source property is the right kind of property.
	if(sourceProperty().containerClass() != subject().dataClass())
		throwException(tr("Modifier was set to operate on '%1', but the selected input is a '%2' property.")
			.arg(subject().dataClass()->pythonName()).arg(sourceProperty().containerClass()->propertyClassDisplayName()));

	// Get modifier inputs.
	const VoxelGrid* voxelGrid = static_object_cast<VoxelGrid>(input.expectLeafObject(subject()));
	voxelGrid->verifyIntegrity();
	OVITO_ASSERT(voxelGrid->domain());
	if(voxelGrid->domain()->is2D())
		throwException(tr("Cannot generate isosurface for a two-dimensional voxel grid. Input must be a 3d grid."));
	const PropertyObject* property = sourceProperty().findInContainer(voxelGrid);
	if(!property)
		throwException(tr("The selected voxel property with the name '%1' does not exist.").arg(sourceProperty().name()));
	if(sourceProperty().vectorComponent() >= (int)property->componentCount())
		throwException(tr("The selected vector component is out of range. The property '%1' contains only %2 values per voxel.").arg(sourceProperty().name()).arg(property->componentCount()));

	for(size_t dim = 0; dim < 3; dim++)
		if(voxelGrid->shape()[dim] <= 1)
			throwException(tr("Cannot generate isosurface for this voxel grid with dimensions %1 x %2 x %3. Must be at least 2 voxels wide in each spatial direction.")
				.arg(voxelGrid->shape()[0]).arg(voxelGrid->shape()[1]).arg(voxelGrid->shape()[2]));

	TimeInterval validityInterval = input.stateValidity();
	FloatType isolevel = isolevelController() ? isolevelController()->getFloatValue(request.time(), validityInterval) : 0;

	// Create engine object. Pass all relevant modifier parameters to the engine as well as the input data.
	return std::make_shared<ComputeIsosurfaceEngine>(validityInterval, voxelGrid->shape(), property->storage(),
			sourceProperty().vectorComponent(), voxelGrid->domain()->data(), isolevel);
}

/******************************************************************************
* Performs the actual analysis. This method is executed in a worker thread.
******************************************************************************/
void CreateIsosurfaceModifier::ComputeIsosurfaceEngine::perform()
{
	setProgressText(tr("Constructing isosurface"));

	if(_mesh.cell().is2D())
		throw Exception(tr("Cannot construct isosurfaces for two-dimensional voxel grids."));
	if(property()->dataType() != PropertyStorage::Float)
		throw Exception(tr("Wrong data type. Can construct isosurface only for floating-point values."));
	if(property()->size() != _gridShape[0] * _gridShape[1] * _gridShape[2])
		throw Exception(tr("Input voxel property has wrong array size, which is incompatible with the grid's dimensions."));

	ConstPropertyAccess<FloatType, true> data(property());
	MarchingCubes mc(_mesh, _gridShape[0], _gridShape[1], _gridShape[2], data.cbegin() + _vectorComponent, property()->componentCount(), false);
	if(!mc.generateIsosurface(_isolevel, *this))
		return;

	// Transform mesh vertices from orthogonal grid space to world space.
	const AffineTransformation tm = cell().matrix() * Matrix3(
		FloatType(1) / (_gridShape[0] - (cell().hasPbc(0)?0:1)), 0, 0,
		0, FloatType(1) / (_gridShape[1] - (cell().hasPbc(1)?0:1)), 0,
		0, 0, FloatType(1) / (_gridShape[2] - (cell().hasPbc(2)?0:1)));
	_mesh.transformVertices(tm);

	// Flip surface orientation if cell matrix is a mirror transformation.
	if(tm.determinant() < 0)
		_mesh.flipFaces();
	if(isCanceled())
		return;

	if(!_mesh.connectOppositeHalfedges())
		throw Exception(tr("Something went wrong. Isosurface mesh is not closed."));
	if(isCanceled())
		return;

	// Determine range of input field values.
	// Only used for informational purposes for the user.
	for(FloatType v : data.componentRange(_vectorComponent)) {
		updateMinMax(v);
	}

	// Compute a histogram of the input field values.
	PropertyAccess<qlonglong> histogramData(histogram());
	FloatType binSize = (maxValue() - minValue()) / histogramData.size();
	int histogramSizeMin1 = histogramData.size() - 1;
	for(FloatType v : data.componentRange(_vectorComponent)) {
		int binIndex = (v - minValue()) / binSize;
		histogramData[std::max(0, std::min(binIndex, histogramSizeMin1))]++;
	}

	// Release data that is no longer needed to reduce memory footprint.
	_property.reset();
}

/******************************************************************************
* Injects the computed results of the engine into the data pipeline.
******************************************************************************/
void CreateIsosurfaceModifier::ComputeIsosurfaceEngine::applyResults(TimePoint time, ModifierApplication* modApp, PipelineFlowState& state)
{
	CreateIsosurfaceModifier* modifier = static_object_cast<CreateIsosurfaceModifier>(modApp->modifier());

	// Look up the input grid.
	if(const VoxelGrid* voxelGrid = dynamic_object_cast<VoxelGrid>(state.expectLeafObject(modifier->subject()))) {
		// Create the output mesh data object.
		SurfaceMesh* meshObj = state.createObject<SurfaceMesh>(QStringLiteral("isosurface"), modApp, tr("Isosurface"));
		mesh().transferTo(meshObj);
		meshObj->setDomain(voxelGrid->domain());
		meshObj->setVisElement(modifier->surfaceMeshVis());
	}

	// Output a data table with the field value histogram.
	DataTable* table = state.createObject<DataTable>(QStringLiteral("isosurface-histogram"), modApp, DataTable::Histogram, modifier->sourceProperty().nameWithComponent(), histogram());
	table->setAxisLabelX(modifier->sourceProperty().nameWithComponent());
	table->setIntervalStart(minValue());
	table->setIntervalEnd(maxValue());

	state.setStatus(PipelineStatus(PipelineStatus::Success, tr("Field value range: [%1, %2]").arg(minValue()).arg(maxValue())));
}

}	// End of namespace
}	// End of namespace
